#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
#-----------------------------------------------------------------------------
# Name:     dsc-noGui.py
# Author:   Cuneyt H. Unlu
# Created:  2013/05/05
# Last mod.:2013/05/21
# Doc:      Conducting thermal analysis using DSC data.
# Prerequisites : Python 2.7, numpy, scipy, matplotlib
#-----------------------------------------------------------------------------
# Changes
# 2013-11-21 Command line interaction was established
#            -p or --pfile argument was ised to process the file
"""
import os
import argparse
import re
import math
import time
import numpy
import scipy
import matplotlib.pyplot as plt
import scipy.interpolate as interp
import scipy.stats as stats


class ols:
    """
    Author: Vincent Nijs (+ ?)

    Email: v-nijs at kellogg.northwestern.edu

    Last Modified: Mon Jan 15 17:56:17 CST 2007

    Dependencies: See import statement at the top of this file

    Doc: Class for multi-variate regression using OLS

    For usage examples of other class methods see the class tests at the bottom of this file. To see the class in action
    simply run this file using 'python ols.py'. This will generate some simulated data and run various analyses. If you have rpy installed
    the same model will also be estimated by R for confirmation.

    Input:
        y = dependent variable
        y_varnm = string with the variable label for y
        x = independent variables, note that a constant is added by default
        x_varnm = string or list of variable labels for the independent variables

    Output:
        There are no values returned by the class. Summary provides printed output.
        All other measures can be accessed as follows:

        Step 1: Create an OLS instance by passing data to the class

            m = ols(y,x,y_varnm = 'y',x_varnm = ['x1','x2','x3','x4'])

        Step 2: Get specific metrics

            To print the coefficients:
                >>> print m.b
            To print the coefficients p-values:
                >>> print m.p

    """

    def __init__(self,y,x,y_varnm = 'y',x_varnm = ''):
        """
        Initializing the ols class.
        """
        self.y = y
        self.x = scipy.c_[scipy.ones(x.shape[0]),x]
        self.y_varnm = y_varnm
        if not isinstance(x_varnm,list):
            self.x_varnm = ['const'] + list(x_varnm)
        else:
            self.x_varnm = ['const'] + x_varnm

        # Estimate model using OLS
        self.estimate()

    def estimate(self):

        # estimating coefficients, and basic stats
        self.inv_xx = numpy.linalg.inv(scipy.dot(self.x.T,self.x))
        xy = scipy.dot(self.x.T,self.y)
        self.b = scipy.dot(self.inv_xx,xy)                    # estimate coefficients

        self.nobs = self.y.shape[0]                     # number of observations
        self.ncoef = self.x.shape[1]                    # number of coef.
        self.df_e = self.nobs - self.ncoef              # degrees of freedom, error
        self.df_r = self.ncoef - 1                      # degrees of freedom, regression

        self.e = self.y - scipy.dot(self.x,self.b)            # residuals
        self.sse = scipy.dot(self.e,self.e)/self.df_e         # SSE
        self.se = numpy.sqrt(numpy.diagonal(self.sse*self.inv_xx))  # coef. standard errors
        self.t = self.b / self.se                       # coef. t-statistics
        self.p = (1-stats.t.cdf(abs(self.t), self.df_e)) * 2    # coef. p-values

        self.R2 = 1 - self.e.var()/self.y.var()         # model R-squared
        self.R2adj = 1-(1-self.R2)*((self.nobs-1)/(self.nobs-self.ncoef))   # adjusted R-square

        self.F = (self.R2/self.df_r) / ((1-self.R2)/self.df_e)  # model F-statistic
        self.Fpv = 1-stats.f.cdf(self.F, self.df_r, self.df_e)  # F-statistic p-value

    def dw(self):
        """
        Calculates the Durbin-Waston statistic
        """
        de = scipy.diff(self.e,1)
        dw = scipy.dot(de,de) / scipy.dot(self.e,self.e);

        return dw

    def omni(self):
        """
        Omnibus test for normality
        """
        return stats.normaltest(self.e)

    def JB(self):
        """
        Calculate residual skewness, kurtosis, and do the JB test for normality
        """

        # Calculate residual skewness and kurtosis
        skew = stats.skew(self.e)
        kurtosis = 3 + stats.kurtosis(self.e)

        # Calculate the Jarque-Bera test for normality
        JB = (self.nobs/6) * (numpy.square(skew) + (1/4)*numpy.square(kurtosis-3))
        JBpv = 1-stats.chi2.cdf(JB,2);

        return JB, JBpv, skew, kurtosis

    def ll(self):
        """
        Calculate model log-likelihood and two information criteria
        """

        # Model log-likelihood, AIC, and BIC criterion values
        ll = -(self.nobs*1/2)*(1+numpy.log(2*numpy.pi)) - (self.nobs/2)*numpy.log(scipy.dot(self.e,self.e)/self.nobs)
        aic = -2*ll/self.nobs + (2*self.ncoef/self.nobs)
        bic = -2*ll/self.nobs + (self.ncoef*numpy.log(self.nobs))/self.nobs

        return ll, aic, bic

    def summary(self):
        """
        Printing model output to screen
        """

        # local time & date
        t = time.localtime()

        # extra stats
        ll, aic, bic = self.ll()
        JB, JBpv, skew, kurtosis = self.JB()
        omni, omnipv = self.omni()

        # printing output to screen
        print '\n=============================================================================='
        print "Dependent Variable: " + self.y_varnm
        print "Method: Least Squares"
        print "Date: ", time.strftime("%a, %d %b %Y",t)
        print "Time: ", time.strftime("%H:%M:%S",t)
        print '# obs:               %5.0f' % self.nobs
        print '# variables:     %5.0f' % self.ncoef
        print '=============================================================================='
        print 'variable     coefficient     std. Error      t-statistic     prob.'
        print '=============================================================================='
        for i in range(len(self.x_varnm)):
            print '''% -5s          % -5.6f     % -5.6f     % -5.6f     % -5.6f''' % tuple([self.x_varnm[i],self.b[i],self.se[i],self.t[i],self.p[i]])
        print '=============================================================================='
        print 'Models stats                         Residual stats'
        print '=============================================================================='
        print 'R-squared            % -5.6f         Durbin-Watson stat  % -5.6f' % tuple([self.R2, self.dw()])
        print 'Adjusted R-squared   % -5.6f         Omnibus stat        % -5.6f' % tuple([self.R2adj, omni])
        print 'F-statistic          % -5.6f         Prob(Omnibus stat)  % -5.6f' % tuple([self.F, omnipv])
        print 'Prob (F-statistic)   % -5.6f			JB stat             % -5.6f' % tuple([self.Fpv, JB])
        print 'Log likelihood       % -5.6f			Prob(JB)            % -5.6f' % tuple([ll, JBpv])
        print 'AIC criterion        % -5.6f         Skew                % -5.6f' % tuple([aic, skew])
        print 'BIC criterion        % -5.6f         Kurtosis            % -5.6f' % tuple([bic, kurtosis])
        print '=============================================================================='


class ImportError(Exception):
    """ Import Error class """
    def __init__(self):
        pass

    def __str__(self):
        return repr("File is not imported. Please check field separator")


class CalcDSC:
    """
    Class for estimating activation energy and kinetic parameters from DSC data.
    """
    def __init__(self,data):
        """ Initialize class with data supplied from user """
        self.importData(data["filepath"]+data["filename"],
                data["field_separator"],
                data["cycle"])

        self.RUN()

    def indices(self, a,  func):
        """ Find index values for values in lists """
        return [i for (i,  val) in enumerate(a) if func(val)]

    def getIndex(self, a, cond):
        """ Get index values for values in lists """
        limitList = self.indices(a, lambda x: x >= cond)
        return limitList[0]

    def importData(self, filename, field_separator,cycle):
        """ Import data.
            filename : used for import data
            field_separator : to determine columns in csv file
            cycle : cycle number of analysis
        """
        infile = open(filename,  "r")

        # read first line and validate if filed_separator is true
        # if it is true go on, else raise ImportError
        line = infile.readline()
        if len(line.split(field_separator)) == 1:
                raise ImportError()
        # check instrument type
        if line.split(field_separator)[0] == 'CLOSED':
            self.instrument = 'TA Q10'
            dismiss = 63
            factor = 1.0
        else:
            print 'Invalid instrument'

        self.cycles = []

        for line in range(dismiss):
            line = infile.readline()
            # find sample weight and heat cycles
            if self.instrument == 'TA Q10':
                if re.search('Sample', line.split(field_separator)[0]):
                    self.sample = str(line.split(field_separator)[1])
                if re.search('Size', line.split(field_separator)[0]):
                    self.sample_weight = float(line.split(field_separator)[1]) / factor
                if re.search('OrgMethod', line.split(field_separator)[0]):
                    if re.search('Mark',line.split(field_separator)[1]):
                        pass
                    else:
                        self.cycles.append(line.split(field_separator)[1])

        time_dsc = []
        temperature = []
        heatflow = []

        time_dsc.append([])
        temperature.append([])
        heatflow.append([])

        k = 0
        while True:
            line = infile.readline()
            if not line:
                break
            tmp = line.split(field_separator)
            if tmp[0] == "-2":
                k += 1
                if k < len(self.cycles):
                    time_dsc.append([])
                    temperature.append([])
                    heatflow.append([])
                pass
            else:
                time_dsc[k].append(float(tmp[0])*60)
                temperature[k].append(float(tmp[1]))
                heatflow[k].append(float(tmp[2])/self.sample_weight)

        self.time = time_dsc[cycle]
        self.temperature = temperature[cycle]
        self.heatflow = heatflow[cycle]

        infile.close()

    def importBaseline(self, filename, field_separator):
        """ if import is True, imports baseline file """
        infile = open(filename,  "r")

        self.baseline_temperature = []
        self.baseline_heatflow = []

        line = infile.readline()

        if len(line.split(field_separator)) == 1:
                raise ImportError()

        if line.split(field_separator)[0] == 'CLOSED':
            self.instrument = 'TA Q10'
            dismiss = 65
        else:
            print 'Invalid instrument'

        for line in range(dismiss):
            line = infile.readline()
            pass

        while True:
            line = infile.readline()
            if not line:
                break
            tmp = line.split(field_separator)

            if str(tmp[0]) == "-2":
                pass
            else:
                self.baseline_temperature.append(float(tmp[1]))
                self.baseline_heatflow.append(float(tmp[2]))
        infile.close()

    def createBaseline(self):
        """ if import is True, creates and subtracts baseline from data """
        # linear interpolation for creating baseline
        self.baseline_heatflow_fit = numpy.interp(self.temperature,self.baseline_temperature,self.baseline_heatflow)

        self.heatflow_fit = []

        dif = \
            self.heatflow[
                        self.decomposition_begin_index] - \
            self.baseline_heatflow_fit[
                        self.decomposition_begin_index]

        for i in range(len(self.temperature)):
            self.heatflow_fit.append(self.heatflow[i] - self.baseline_heatflow_fit[i] - dif)

    def integrateHeatflow(self, idx_1, idx_2):
        """ integrates heatflow against time for a given range """
        return numpy.trapz(self.heatflow_data[idx_1:idx_2], self.time_data[idx_1:idx_2])

    def findPercent(self,percent,data):
        """ determination of signals having given percent of main signal """
        for i in range(len(data)):
            proportion = data[i] / self.max_heatflow
            if proportion > percent - 0.0001 and proportion < percent + 0.01:
                value = data[i]
                value_index = i
                break
        return value, value_index

    def createAxes(self, first, second, step):
        """ creation of axes for calculation """
        # alpha (a)
        self.alpha = []
        # ln (da/dt)
        self.ln_d_alpha_d_t = []
        # 1/T
        self.inverse_temperature = []
        # ln (a)
        self.ln_alpha = []
        for i in range(first, second, step):
            # 1 - alpha = delta_h / delta_h_total
            self.alpha.append(self.integrateHeatflow(i, -1) / self.delta_h_total)
            # d_alpha/d_t = (d_H/d_t)/delta_h_total
            self.ln_d_alpha_d_t.append(math.log(self.heatflow_data[i] / self.delta_h_total))
            # 1/T
            self.inverse_temperature.append(1 / (self.temperature[self.decomposition_begin_index + i] + 273.15))
            # ln(1-alpha)
            self.ln_alpha.append(math.log(self.integrateHeatflow(i, -1) / self.delta_h_total))

    def fitLine(self,x,y):
        """ linear fit function """
        tmp = numpy.polyfit (x, y, 1)
        p = numpy.poly1d(tmp)
        return tmp, p

    def plots(self):
        """ Plots of data """
        # main figure
        figure = plt.figure(1)
        # heatflow / temperature
        plot_dsc = figure.add_subplot(221)
        # derivative heatflow/temperature
        plot_derivative = figure.add_subplot(223)
        # calculation range
        plot_calc = figure.add_subplot(222)
        # ln k / 1/T
        plot_ln_k = figure.add_subplot(224)

        plot_dsc.plot(self.temperature,
                      self.y_data,
                      'k-',
                      linewidth=2)

        plot_ln_k.plot(self.inverse_temperature,
                       self.ln_k,
                       'ko',
                       self.inverse_temperature,
                       self.ln_k_fit,
                       'r-',
                       linewidth=2)

        plot_derivative.plot(self.temperature[1:],
                             self.derivative_dsc,
                             'k-',
                             linewidth=2)
        plot_calc.plot(self.time_data,
                       self.heatflow_data,
                       'k-',
                       self.time_data[self.peak_height_1_index:self.peak_height_2_index],
                       self.heatflow_data[self.peak_height_1_index:self.peak_height_2_index],
                       'r-',
                       linewidth=2)

        # Plot formatting
        plot_dsc.set_title('Heatflow')
        plot_dsc.set_ylabel('Heatflow (mW/g)')
        plot_dsc.grid(True)

        plot_ln_k.set_xlabel('1000/T')
        plot_ln_k.set_ylabel('ln(k) (1/s)')
        plot_ln_k.grid(True)
        plot_ln_k.xaxis.get_major_formatter().set_powerlimits((-3,0))
        plot_ln_k.xaxis.get_major_formatter().set_scientific(True)

        plot_derivative.set_xlabel('Temperature (C)')
        plot_derivative.set_ylabel('Derivative Heatflow (mW/g C)')
#        plot_derivative.set_ylim((-0.1,0.1)) % for y limits to display
        plot_derivative.grid(True)

        plot_calc.set_title('Calculation range')
        plot_calc.set_xlabel('Time (s)')
        plot_calc.set_ylabel('Heatflow (mW/g)')
        plot_calc.grid(True)

        plt.tight_layout()

        if data["plots"]:
            plt.show()

        if data["report"]:
            plt.savefig(self.sample+'.png',  dpi=100,  bbox_inches='tight')

    def displayResults(self):
        """ prints results to the console """
        print '-'*60
        print "Date of Analysis:\t %s" % self.timestamp
        print "Sample Name: \t %s" % self.sample
        print "Sample \n \tPath: \t %s\n \tFilename: \t %s" % (data["filepath"],
                                                               data["filename"])
        if data["baseline"]:
            print "Baseline \n \tPath: \t %s\n \tFilename: \t %s" % (data["baseline_path"],
                                                                     data["baseline_filename"])
        print '-'*60
        print "Decomp. Limits \tBegin: %4.2f \t End: %4.2f " % (data["decomposition_begin"],data["decomposition_end"])
        print '-'*60
        print "Decomposition heat \t %4.2f J/g" % self.delta_h_total
        print "Max. Decomposition \t %4.2f C" % self.temperature[self.max_heatflow_index]
        print "%d%% peak position \t %3.1f s \t %4.2f C" % (self.calculation_peak_height_1*100,
                                                         self.time_data[self.peak_height_1_index],
                                                         self.temperature[self.decomposition_begin_index+self.peak_height_1_index])
        print "%d%% peak position \t %3.1f s \t %4.2f C" % (self.calculation_peak_height_2*100,
                                                         self.time_data[self.peak_height_2_index],
                                                         self.temperature[self.decomposition_begin_index+self.peak_height_2_index])
        print "Sampling data for every %d value " % self.step
        print '-'*60
        print "Regression results"
        print "ln Z \t %4.2f 1/s" % self.ln_z
        print "Order \t %3.1f" % self.order
        print "Ac.En. \t %4.2f kJ/mol" % self.activation_energy
        print "R sq. \t %6.5f" % self.m.R2
        print '-'*60
        print "1000/T \t T (K) \t t (C) \t ln(k) \t k*1000 (1/s)"
        print '-'*60
        for i in range(len(self.inverse_temperature)):
            print "%4.3f \t %4.1f \t %4.1f \t %4.2f \t %4.3f" % (1000 * self.inverse_temperature[i],
                                      1/self.inverse_temperature[i],
                                      1/self.inverse_temperature[i]-273.15,
                                      self.ln_k[i],
                                      1000 * math.exp(self.ln_k[i]))
        print '-'*60

    def createReport(self):
        """ Creates an HTML report """
        report_file = open(self.sample+".html",'w')
        report_file.write(
        """
            <!DOCTYPE html>
            <style type="text/css">
            tr:nth-of-type(odd) {
              background-color:#ccc;
            }
            table.table_class {
            	border-width: 0px;
            	border-spacing: 0px;
            	border-style: none;
            	border-collapse: separate;
            }
            table.table_class th {
            	border-width: 0px;
            	padding: 0.2em 2ex 0.2em 2ex;
            	border-style: none;
            }
            table.table_class td {
            	border-width: 0px;
            	padding: 0.2em 2ex 0.2em 2ex;
            	border-style: none;
            	border-color: gray;
            }
            </style>

            <html>
            <head>
            <meta charset="utf-8" />
            <title>Report of the activation energy calculation</title>
            </head>
            <body>
        """)
        report_file.write(""" <table> """)
        report_file.write(
        """
        <h1><b>%s <!-- sample name --></b></h1></br>
        <tr>
            <td><h4><em>Time of Analysis</em></h4> </td>
            <td>%s <!-- date and time --></td>
        </tr>
        <tr>
            <td><h4>File path</h4></td>
            <td>%s <!-- file path --></td>
        </tr>
        <tr>
            <td><h4>File name</h4></td>
            <td>%s <!-- file name --></td>
        </tr>
        """
        %
        (self.sample,
         self.timestamp,
         data["filepath"],
         data["filename"]
         )
        )
        if data["baseline"]:
            report_file.write(
            """
            <tr>
                <td><h4>Baseline file path</h4></td>
                <td>%s <!-- file name --></td>
            </tr>
            <tr>
                <td><h4>Baseline file name</h4></td>
                <td>%s <!-- file name --></td>
            </tr>
            """
            %
            (data["baseline_path"],
             data["baseline_filename"]
             )
            )
        report_file.write("""
            <tr>
                <td><h4>Instrument </h4></td>
                <td>%s <!-- instrument --></td>
            </tr>
            """
            %
            (
            self.instrument
            )
        )
        report_file.write("""</table> """)
        report_file.write("""<table> """)
        report_file.write(
        """
            <table class="table_class">
            <tr>
            	<th>Temperature limits (<sup>o</sup>C)</th>
            	<td><b>Begin</b></td>
            	<td><b>End</b></th>
            </tr>
            <tr>
            	<td>Decompostion</td> <!-- Decomposition limits -->
            	<td> %5.2f </td>
                <td> %5.2f </td>
            </tr>
            <tr>
            	<th>Max. degradation</th>
            	<td>%5.2f</td>
            	<td></td>
            </tr>
            <tr>
            	<th>Decomposition Heat (&Delta;Q) </th>
            	<td>%5.2f</td>
            	<td> J/g </td>
            </tr>
        """
        %
        (
            self.temperature[
                    self.decomposition_begin_index
                    ],
            self.temperature[
                    self.decomposition_end_index
                    ],
            self.temperature[self.max_heatflow_index],
            self.delta_h_total
            )

        )
        report_file.write(
        """
        <tr>
            <td>%d%% peak position </td>
            <td> %3.1f s </td>
            <td> %4.2f C </td>
        </tr>
        """ % (
            self.calculation_peak_height_1*100,
            self.time_data[self.peak_height_1_index],
            self.temperature[self.decomposition_begin_index+self.peak_height_1_index]
            )
        )
        report_file.write(
        """
        <tr>
            <td>%d%% peak position </td>
            <td> %3.1f s </td>
            <td> %4.2f C </td>
        </tr>
        """ % (
            self.calculation_peak_height_2*100,
            self.time_data[self.peak_height_2_index],
            self.temperature[self.decomposition_begin_index+self.peak_height_2_index]
            )
        )
        report_file.write(
        """
        <tr>
            <td> Sampling data for every </td>
            <td> %d </td>
            <td> value </td>
        </tr>
        """ % self.step
        )
        report_file.write("""</table> """)
        report_file.write( """ <h3> Regression results </h3>""")
        report_file.write("""<table> """)
        report_file.write(
        """
        <tr>
            <td width="50"> ln Z </td>
            <td width="60"> %4.2f </td>
            <td width ="175"> Arrhenius factor (1/s) </td>
            <td width="50"> Order </td>
            <td width="60"> %3.1f </td>
            <td width="175"> Reaction order </td>
        </tr>
        <tr>
            <td> E<sub>A</sub> </td>
            <td> %4.2f </td>
            <td> Activation Energy (kJ/mol) </td>
            <td> R </td>
            <td> %.5f </td>
            <td> Correlation coefficient </td>
        </tr>
        """
        %
            (
             self.ln_z,
             self.order,
             self.activation_energy,
             self.m.R2
            )
        )
        report_file.write("""</table> """)
        report_file.write("""<table> """)
        report_file.write(
        """
            <tr>
            <th> Inv. temp. (1000/K) </th>
            <th> Abs. temp. (K) </th>
            <th> Temp. (C) </th>
            <th> ln (k) (1/s) </th>
            <th> 1000 x k (1/s) </th>
            </tr>
        """)
        for i in range(len(self.inverse_temperature)):
            report_file.write(
            """
            <tr>
            <td><center> %4.3f <center></center> </td>
            <td><center> %4.1f <center></td>
            <td><center> %4.1f <center></td>
            <td><center> %4.3f <center></td>
            <td><center> %4.3f <center></td>
            </tr>
            """
            %
            (1000 * self.inverse_temperature[i],
             1/self.inverse_temperature[i],
             1/self.inverse_temperature[i]-273.15,
             self.ln_k[i],
             1000 * math.exp(self.ln_k[i])
            )
        )
        report_file.write("""</table> """)
        report_file.write(
            """
            <Img src="%s.png" width="600px" alt="" />
            </body>
            </html>
            """
            %
            self.sample
            )

        report_file.close()

    def exportData(self):
        """ exports 1/T and ln k data """
        export_file = open(self.sample+".XY","w")
        export_file.write("# Inverse temperature, ln k,  data\n")
        export_file.write("# Sample %s\n" % self.sample)
        export_file.write("# Date %s \n" % self.timestamp)
        for i in range(len(self.inverse_temperature)):
           export_file.write("%6.5f \t %6.5f"
           %(
            self.inverse_temperature[i],
            self.ln_k[i]
            )
        )
        export_file.close()

    def RUN(self):
        """ run calculation """
        self.timestamp = time.strftime("%d-%b-%Y %H:%M:%S", time.localtime())

        self.decomposition_begin = data["decomposition_begin"]
        self.decomposition_end = data["decomposition_end"]

        self.calculation_peak_height_1 = data["calculation_peak_height_1"]
        self.calculation_peak_height_2 = data["calculation_peak_height_2"]

        self.step = data["step"]
        # indices for starting of decomopisition
        self.decomposition_begin_index = self.getIndex(
                                                    self.temperature,
                                                    self.decomposition_begin)
        self.decomposition_end_index = self.getIndex(
                                                    self.temperature,
                                                    self.decomposition_end)

        if data["baseline"]:
            self.importBaseline(data["baseline_path"]+data["baseline_filename"],
                                data["field_separator"])
            self.createBaseline()
            self.y_data = self.heatflow_fit
        else:
            self.y_data = self.heatflow

        self.max_heatflow = max(self.y_data)

        self.max_heatflow_index = self.getIndex(
                                                self.y_data,
                                                self.max_heatflow)

        self.time_data = [x - self.time[self.decomposition_begin_index] for x in self.time[self.decomposition_begin_index:self.decomposition_end_index]]
        self.heatflow_data = self.y_data[self.decomposition_begin_index:self.decomposition_end_index]

        self.peak_height_1, self.peak_height_1_index = self.findPercent(self.calculation_peak_height_1, self.heatflow_data)
        self.peak_height_2, self.peak_height_2_index = self.findPercent(self.calculation_peak_height_2, self.heatflow_data)

        self.delta_h_total = self.integrateHeatflow(0, -1)

        self.createAxes(self.peak_height_1_index, self.peak_height_2_index, self.step)
        # multiple linear regression
        z = self.ln_d_alpha_d_t
        x = self.ln_alpha
        y = self.inverse_temperature

        X = numpy.array([x,y]).T
        # multiple linear regression function
        self.m = ols(numpy.transpose(z),X,y_varnm = 'z',x_varnm = ['x', 'y'])
        # m is coefficents vector
        self.ln_z = self.m.b[0]
        self.order = self.m.b[1]
        self.activation_energy = self.m.b[2]/8.314/1000

        self.ln_k = [self.ln_z + self.m.b[2] * i for i in y ]

        tmp, p = self.fitLine(self.inverse_temperature, self.ln_k)

        self.ln_k_fit = [p(z) for z in self.inverse_temperature]

        self.derivative_dsc = list(numpy.diff(self.y_data)/numpy.diff(self.temperature))

def runCalculation():
    """ Calculation function of the module. Takes no variable """
    global result
    # Don't change anything here
    # if you don't know what you are doing exactly,
    # otherwise script will not work correctly!
    data.update()
    # create Calc instance
    result = CalcDSC(data)
    result.plots()

    if data["screen"]:
        result.displayResults()

    if data["report"]:
        result.createReport()

    if data["export"]:
        result.exportData()

def fileData():
    """ user input part """
    global data
    data = {}
    # sample parameter file
    data["filepath"] = os.getcwd() + os.sep

    parser = argparse.ArgumentParser()

    parser.add_argument('--pfile'   , action='store', dest='pfile')
    parser.add_argument('-p'   , action='store', dest='pfile')

    args = parser.parse_args()

    if args.pfile is None:
        # Argument `--C` not given on command line
        args.pfile = "pva.txt"
    execfile(args.pfile)


#   ## File name and other variables
#   ## i.e.
#   ## data["filepath"] = "/home/user/data/"
#   ## or
#   ## data["filepath"] = "C:\\data\\"
#   ## by default this is the folder where script runs
#   data["filepath"] = os.getcwd() + os.sep
#   ## Change data file name
#   ## sample files
#   ## PAN homopolymer recorded on TA Q10
#   data["filename"] = "DSC_sample_PAN.csv"
#   ## Similar behavior is also present for baseline file
#   ## if there are any baseline file change Boolean value
#   ## and give the name of the file for baseline correction
#   data["baseline"] = True
#   if data["baseline"]:
#       data["baseline_path"] = os.getcwd() + os.sep
#       data["baseline_filename"] = "DSC_sample_baseline.csv"
#
#    ##  Select data file separator,  for example \t for tab,  etc.
#    data["field_separator"] = ","
#    ## calculation parameters
#    ## This is input section of calculation variables
#    ## start - end points of decomposition
#   data["decomposition_begin"] = 295
#   data["decomposition_end"] = 326
#   ## peak heights for calculation
#   data["calculation_peak_height_1"] = 0.15
#   data["calculation_peak_height_2"] = 0.80
#   ## calculation sampling step for calculation
#   data["step"] = 10
#   ## Heat cycle to be used or calculation
#   ## default is last cycle
#   data["cycle"] = -1
#
#   ## This section includes Boolean (True/False) values
#   ## Do you want to see results on screen ?
#   data["screen"] = True
#   ## Do you want to see plots on screen ?
#   data["plots"] = True
#   ## Do you want a report as html file ?
#   data["report"] = False
#    ## Do you want to export data as ASCII files?
#    data["export"] = False


if __name__ == "__main__":
    # retrieve filename and other variables
    fileData()
    # make the calculation
    runCalculation()

# vim: tabstop= 4 expandtab shiftwidth= 4 softtabstop= 4
